setwd("___YOUR FILE PATH HERE___") 


rm(list=ls())
set.seed(12435)
library(foreign)
library(interplot)
library(ltm)
library(stargazer)
library(tidyverse)
library(lmtest)
library(multiwayvcov)
library(plm)
library(gridExtra)
library(emmeans)
library(broom)
library(dotwhisker)
require(lattice)
library(reshape2)
library(sjPlot)
library(effects)
library(ape)
library(dplyr)

devtools::install_github('gibbonscharlie/bfe')
library(bfe)


# read in the data
elect_data<-read.csv('elections_2017.csv', stringsAsFactors = F)

# need to delete states with no internal variation in Distance for the reweighting analysis
elect_data_bfe <- subset(elect_data, NAME_1!="Berlin" & NAME_1!="Hamburg")


#### BASELINE ANALYSES, WITHOUT THEN WITH FEs
mod1a <- lm(AfDshare ~ distance2, elect_data)
mod2a <- lm(AfDNPDshare ~ distance2, elect_data)

mod1b <- lm(AfDshare ~ distance2 + I(NAME_1), elect_data)
mod2b <- lm(AfDNPDshare ~ distance2 + I(NAME_1), elect_data)

mod1c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden, elect_data)
mod2c <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden, elect_data)

mod1d <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + I(NAME_1), elect_data)
mod2d <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + I(NAME_1), elect_data)

tab_model(mod1a, mod1b, mod1c, mod1d,digits=3,
          show.se=TRUE, show.ci=FALSE, collapse.se=TRUE, show.p=FALSE, p.style="asterisk", p.threshold = c(0.05, 0.01))

tab_model(mod2a, mod2b, mod2c, mod2d,digits=3,
          show.se=TRUE, show.ci=FALSE, collapse.se=TRUE, show.p=FALSE, p.style="asterisk", p.threshold = c(0.05, 0.01))


#########  HAUSMAN TESTS
form <- AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden
wi<- plm(form, data = elect_data, model = "within", index="NAME_1")
re <- plm(form, data = elect_data, model = "random", index="NAME_1")
pooled <- plm(form, data = elect_data, model = "pooling", index="NAME_1")
re.v.pooled <- round(phtest(re, pooled)$p.value, digits=3)
fe.v.pooled <- round(phtest(wi, pooled)$p.value, digits=3)
fe.v.re <- round(phtest(wi, re)$p.value, digits=3)

form <- AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden
wi2<- plm(form, data = elect_data, model = "within", index="NAME_1")
re2 <- plm(form, data = elect_data, model = "random", index="NAME_1")
pooled2 <- plm(form, data = elect_data, model = "pooling", index="NAME_1")
re.v.pooled2 <- round(phtest(re2, pooled2)$p.value, digits=3)
fe.v.pooled2 <- round(phtest(wi2, pooled2)$p.value, digits=3)
fe.v.re2 <- round(phtest(wi2, re2)$p.value, digits=3)


stargazer(pooled,re,wi, pooled2,re2,wi2, 
          out="table6.txt", no.space=T, keep="distance2",
          star.cutoffs = c(0.05, 0.01, NA), keep.stat = c("n","adj.rsq"),
          add.lines=list(c("Method","Pooled","RE","FE","Pooled","RE","FE"),
                         c("RE v Pooled",re.v.pooled,"","",re.v.pooled2,"",""),
                         c("FE v Pooled",fe.v.pooled,"","",fe.v.pooled2,"",""),
                         c("FE v RE",fe.v.re,"","",fe.v.re2,"","")))





###### IWE AND RWE ANALYSES; CAN'T ADD TO STARGAZER BUT CAN STORE TO ADD MANUALLY

mod1e <- EstimateIWE(y="AfDshare", treatment="distance2", group="NAME_1", 
                     controls=c("prop_juden","unemp33","pop25","vshare33"),
                     data=elect_data_bfe, subset=NULL, cluster.var=NULL,is.robust=TRUE)
mod1f <- EstimateRWE(y="AfDshare", treatment="distance2", group="NAME_1", 
                     controls=c("prop_juden","unemp33","pop25","vshare33"),
                     data=elect_data_bfe, subset=NULL, cluster.var=NULL,is.robust=TRUE)
m1.iwe <- rbind(mod1e$swe.est, sqrt(mod1e$swe.var), mod1e$swe.est/sqrt(mod1e$swe.var))
m1.rwe <- rbind(mod1f$swe.est, sqrt(mod1f$swe.var), mod1f$swe.est/sqrt(mod1f$swe.var))

res.m1 <- cbind(m1.iwe, m1.rwe)
colnames(res.m1) <- c("IWE","RWE")

mod2e <- EstimateIWE(y="AfDNPDshare", treatment="distance2", group="NAME_1", 
                     controls=c("prop_juden","unemp33","pop25","vshare33"),
                     data=elect_data_bfe, subset=NULL, cluster.var=NULL,is.robust=TRUE)
mod2f <- EstimateRWE(y="AfDNPDshare", treatment="distance2", group="NAME_1", 
                     controls=c("prop_juden","unemp33","pop25","vshare33"),
                     data=elect_data_bfe, subset=NULL, cluster.var=NULL,is.robust=TRUE)
m2.iwe <- rbind(mod2e$swe.est, sqrt(mod2e$swe.var), mod2e$swe.est/sqrt(mod2e$swe.var))
m2.rwe <- rbind(mod2f$swe.est, sqrt(mod2f$swe.var), mod2f$swe.est/sqrt(mod2f$swe.var))

res.m2 <- cbind(m2.iwe, m2.rwe)
colnames(res.m2) <- c("IWE","RWE")


write.table(round(cbind(res.m1,res.m2),3), file="table6_columns4and5.txt", row.names=FALSE)



#### TWOWAY FES
mod1k <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + I(closest_camp2) + I(NAME_1), elect_data)
mod2k <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + I(closest_camp2) + I(NAME_1), elect_data)

#### TWOWAY RES
mod1l <- lmer(AfDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                          ( 1 | NAME_1) + (1 | closest_camp2),elect_data,
                        control = lmerControl(optimizer ="Nelder_Mead"))
mod2l <- lmer(AfDNPDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                          ( 1 | NAME_1) + (1 | closest_camp2),elect_data,
                        control = lmerControl(optimizer ="Nelder_Mead"))

### 3 way random effects by camp, Land, Kreis

mod1m <- lmer(AfDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                ( 1 | NAME_1) + (1 | NAME_2) + (1 | closest_camp2),elect_data,
              control = lmerControl(optimizer ="Nelder_Mead"))
mod2m <- lmer(AfDNPDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                ( 1 | NAME_1) + (1 | NAME_2) + (1 | closest_camp2),elect_data,
              control = lmerControl(optimizer ="Nelder_Mead"))


#### TWOWAY RANDOM COEFFICIENTS
mod1n <- lmer(AfDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                ( 1  + distance2 | NAME_1) + (1  + distance2 | closest_camp2),elect_data,
              control = lmerControl(optimizer ="Nelder_Mead"))
mod2n <- lmer(AfDNPDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                ( 1 + distance2 | NAME_1) + (1 + distance2 | closest_camp2),elect_data,
              control = lmerControl(optimizer ="Nelder_Mead"))

### 3 way random coefficients!!!!

mod1o <- lmer(AfDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                (1 + distance2 | NAME_1) + 
                (1 + distance2 | NAME_2) + 
                (1 + distance2 | closest_camp2),elect_data,
              control = lmerControl(optimizer ="Nelder_Mead"))
mod2o <- lmer(AfDNPDshare~distance2 + vshare33 + unemp33 + pop25 + prop_juden + 
                (1 + distance2 | NAME_1) + 
                (1 + distance2 | NAME_2) + 
                (1 + distance2 | closest_camp2),elect_data,
              control = lmerControl(optimizer ="Nelder_Mead"))


tab_model(mod1k,mod1l,mod1m,mod1n,mod1o,mod2k,mod2l,mod2m,mod2n,mod2o,digits=3, terms=c("distance2"), 
          show.se=TRUE, show.ci=FALSE, collapse.se=TRUE, show.p=FALSE, p.style="asterisk", p.threshold = c(0.05, 0.01))
stargazer(mod1k,mod1l,mod1m,mod1n,mod1o,mod2k,mod2l,mod2m,mod2n,mod2o, dep.var.labels.include = FALSE,type="text",
          out="table7.txt", no.space=T, 
          keep = c("distance"),
          star.cutoffs = c(0.05, 0.01, NA),  keep.stat = c("n","adj.rsq"),
          add.lines=list(c("Method","Fixed Effects","Random Intercept","Random Coefficient","Fixed Effects","Random Intercept","Random Coefficient")), 
          covariate.labels = c('Distance to camp'))

anova(mod1l,mod1m,mod1n,mod1o)
anova(mod2l,mod2m,mod2n,mod2o)



